Reconstructing past migratory behaviour of reindeer (Rangifer tarandus): Insights from geometric morphometric analysis of proximal phalanx morphology from extant caribou populations

Reindeer mobility patterns vary widely in modern ecosystems, notably between more open or more wooded environments. This renders the reconstruction of past reindeer mobility patterns challenging, while being at the same time key if archaeologists want to better understand the role that reindeer herds played in the subsistence and territorial organisation of Prehistoric hunter-gatherer societies. Studying the morphology associated with different habitats and mobility patterns can be a useful method for understanding past reindeer behaviour. To access paleoecological information, the relationship between locomotor anatomy and substrate type can be explored in modern animals and transposed to the past, as previous studies have proven that an animal´s environment affects bone morphology. In this study, 3D Geometric Morphometrics are used to explore the impact of extant reindeer habitat type and mobility pattern on phalanx morphology. Results obtained reflects on the potential archaeological application of such an approach for paleoecological reconstructions. Size and shape of phalanx vary significantly, yet complex to interpret in light of interplaying factors such as subspecies, sexual dimorphism and the influence of migration costs, snow cover and substrate type. If direct application to the archaeological record remains preliminary, this first study permits to highlight promising avenues for future research.


Introduction
During the Late Pleistocene, the reindeer (Rangifer tarandus) was one of the key prey species for many Palaeolithic hunter-gatherer groups [e.g. [1][2][3]. In Western Europe, several archaeological sequences documented hunting episodes focused on reindeer herds and whose studies and identification in the archaeological record (particularly to distinguish the different stages of its domestication). They have provided an important understanding of this subject, which is key in reconstructing the mobility and lifestyles of northern peoples [68,[70][71][72][73][74][75][76][77]. In parallel, the increased development of geometric morphometrics (GMM) provides a powerful set of techniques for distinguishing subtle morphological differences among members of the same species. It was successfully applied to the study of ecomorphological patterns in several instances [67,[75][76][77][78][79][80].
Building from these previous works, we here extend the application of 3D GMM to phalanges and to a larger sample of wild sedentary and (latitudinal or altitudinal) migratory populations, in order to explore the potential link between proximal phalanges morphology and reindeer mobility and habitat. The choice of the proximal phalanx is due to the fact that it is a often element in the archaeological record of the potential sites to be studied.
Thus, the main purpose of this work is to determine if there are size and shape differences in proximal phalanx morphology between subspecies (Rangifer tarandus caribou, Rangifer tarandus granti, Rangifer tarandus groenlandicus, Rangifer tarandus pearyi) that have different mobility patterns (altitudinal, longitudinal, sedentary) and live in different habitats (mountain, boreal forest and tundra). Moreover, since Rangifer tarandus caribou subspecies lives in mountainous areas, boreal forests, tundra, and undertakes altitudinal and longitudinal migrations or has a sedentary behaviour [43,81], this subspecies (see Sampling section for further details) will also be studied separately.
Although morphological studies of reindeer often consider sex a key variable [75,76,82] or even other ungulate osteometric studies [83,84], the presence of sexual dimorphism will be just visually explored in size (distribution of the log-centroid sizes) in the subspecies Rangifer tarandus caribou for which there is a limited but available sample of males and females.
Finally, as previous studies have shown that allometry is significant [77], this study will explore the effect of size on the shape of proximal phalanges of both forelimbs and hindlimbs.

Sampling extant caribou populations
Rangifer tarandus has a wide circumpolar distribution. For this study, we acquired proximal phalanges of individuals belonging to the four subspecies of North American caribou (see Fig 1) (mainly from Canada). These subspecies occupy different habitat settings (boreal, montane and arctic environments), sometimes overlapping, and adopt different mobility strategies. Across this large area, caribou exhibit a wide variation in ecology, genetics, behaviour and morphology. Biologists generally distinguish two main mobility patterns for caribou: migratory and sedentary [42]. However, in mountainous settings some caribou populations perform these migratory movements altitudinally [85], calling it altitudinal or elevational migration.
The subspecies selected for this study are the following: 1) R. t. caribou (which includes Eastern caribou/Migratory woodland caribou, mountain caribou and boreal caribou/woodland caribou) inhabits different environments and performs short, medium and long-distance movements. In Quebec, Woodland caribou living south of 55˚N are forest-dwellers and seasonally move relatively short distances compared to their tundra counterparts [86,87], between 50 and 150 kilometres depending on the herd [49,48,88,89]. Eastern or Migratory Woodland caribou group two big herds, George River and Leaf River, that, despite belonging to Woodland caribou, behave like Barren-Ground or tundra caribou, performing migration on long distances ranging between 1120-1770 kilometres on average per year [90], from open tundra settings to boreal forest habitats [91]. Mountain caribou inhabits mountain and boreal settings. Some herds perform altitudinal migrations, with overall distances between 140-240 kilometres only in spring migration, depending upon the location of their calving areas, while others remain relatively sedentary, travelling shorter distances under 100 kilometres [81,92].

PLOS ONE
subspecies [69]. The first one is mainly found in Alaska and its adjacent Canadian territories, the second one in Canada only (Nunavut, Northwest Territories). These tundra caribou include some of the largest herds in North America, numbering hundreds of thousands of individuals that make very long-distance migrations, ranging between 800 and 5055 kilometres (annual migration) depending on the herds and the year [81,[100][101][102][103]. Their habitats are mainly prostrate dwarf shrub tundra and upright shrub tundra [94].
For this research we compiled a sample of 34 fore-and 44 hindlimb proximal (first) phalanges belonging to the four subspecies of caribou described above (Tables 1 and 2). All of them were adults and specimens with bone pathologies were excluded from the study. Samples were obtained from biologists in the Ministry of Forests, Lands, Natural Resource Operations and Rural Development, British Columbia, who provided us with mountain caribou samples and related information about herds (type of herd, sex if known and location). Peary caribou, Barren-Ground and Grant's caribou specimens were collected from the Canadian Museum of Nature (Ottawa), which owns one of the largest collections of caribou in the country. Eastern migratory and short distance woodland caribou were obtained with the help of Biologists from the Forest and Wildlife branch of the Québec Government, the Prehistory and bioarchaeology lab from Laval University and the archaeozoology laboratories of the University of Montréal.
To facilitate data analysis, each individual was given a unique number (from DPpha 1 to DPpha X, with DP standing for "DeerPal", the acronym of the project). A detailed table listing the habitat and mobility pattern of each phalanx, the number of our unique ID, as well as population information is provided in the S1 Table. The fact that sedentary populations are classified as "endangered" by Canadian Government makes it more difficult to obtain samples from these herds (since there are very few specimens in museums, the only way to obtain them was through samples sent by biologists when a caribou was killed in an accident or found dead) resulting in uneven samples. Unfortunately, sex was unknown for almost all individuals (n = 14 forelimb phalanges; n = 18 hindlimb phalanges).
Groups were considered according to Subspecies (Rangifer tarandus caribou, Rangifer tarandus granti, Rangifer tarandus groenlandicus, Rangifer tarandus pearyi), Mobility (Altitudinal for populations migrating vertically, Planarly for populations migrating mostly on a latitudinal or longitudinal plane, and Sedentary) and Habitat (Mountain, Boreal Forest and Tundra) (see Table 1). In order to analyse R.t. caribou subspecies, specimens were grouped according to habitat and mobility pattern (see Table 3).

Data acquisition and 3D geometric morphometrics
Phalanges were individually scanned using a NextEngine 3D laser Scanner (at the Anthropology Department of the Université de Montréal, Canada), in order to obtain a highly accurate Geometric morphometrics (GMM) is a quantitative approach which allows the comparison of bone shapes and the visualization of significant morphological differences between groups of specimens while retaining the element of shape information related to size [104,105].
Because of the uncertainty in terms of the anatomical position of the curve and surface semilandmarks, sliding step was carried out in curves along their tangent vectors and surface points within their tangent planes in order to minimize the bending energy of the landmark configurations and the template during the Thin-Plate Spline (TPSE) interpolation [106,107]. After that, the resulting landmark configurations were rotated, translated and scaled using the Generalized Procrustes Analysis (GPA) [108,109]. This analysis eliminates the effect of position and orientation of the configurations, and limits the effect of size [110][111][112], yielding Procrustes shape coordinates that were subsequently analysed.
Considering that proximal phalanges are often found fragmented in Palaeolithic sites, we focused our work on proximal epiphyses, in order to allow future applications to a larger part of the archaeological record. The round shape of articular surfaces of these phalanges make it difficult to define true anatomical landmarks, and thus semilandmarks were preferred to define the complex morphology of the joint [106]. We placed landmarks and semilandmarks on areas not affected by entheseal changes and paleopathologies [70,72]. Two landmark configurationtemplates were made with a total of 4 landmarks (LM), 43 curve semilandmarks (CSLM) and 25 surface semilandmarks (SSLM) for forelimb phalanges and 4 LM, 45 CSLM and 25 SSLM for hindlimb phalanges, that were digitized to get the form of the proximal epiphyses (Fig 2, Tables 4 and 5). Digitalization was carried out using ViewBox 4.1.0.12 (dHAL software, Kifissia, Greece).

Statistical analyses: Size and shape differences
First, the four subspecies were analysed altogether. Then, since R. t. caribou subspecies constituted a greater percentage (70.5%-forelimb phalanges and 81.8%-hindlimb phalanges) of the sample, the same analyses were only conducted on individuals from this subspecies.
Thus, based on log-transformed centroid sizes for forelimb and hindlimb phalanges, the size differences were assessed using Kruskal-Wallis tests (threshold set at p< 0.05) for the specimens pooled by "Subspecies", "Habitat", "Mobility", and the interactions "Subspecies https://doi.org/10.1371/journal.pone.0285487.g002 +Habitat", "Subspecies+Mobility" and "Mobility+Habitat". Then, multiple Wilcoxon-rank tests were used to compare these groups according to these categories. P values were corrected using the "Benjamini-Hochberg" method to decrease the false discovery rate [113].
A Principal Component Analysis (PCA) based on Procrustes Coordinates was conducted to observe shape variation in the morphospace. Using a Thin-Plate Spline (TPS) interpolation function [108] we generated a 3D digital mesh for each articular surface respectively that was wrapped toward the Procrustes grand mean to understand variations along the principal axes. From the surface of the Procrustes mean configurations [114], shapes at the extremes of the axes of the principal components were visualized and magnified by a scale factor of 0.1 (in order to improve the visualization of the morphological differences), using the function "tps3d" in the "Morpho" package [115]. A multivariate analysis of variance (MANOVA) on PC scores was used to estimate shape differences between these groups. Additionally, in order to test the assignment accuracy for each category, a Canonical Variate Analysis (CVA) was computed on PC scores (data dimensionality reduction) by calculating the cross-validated correct classification and the Kappa metric (see S1 File). We avoided overfitting [116] by never exceeding 95% of the variance cumulative and by never using more PCs than the number of specimens from the smallest group or category.
As sex variable was known in some specimens from R. t. caribou subspecies, sexual size dimorphism was only explored in this subspecies through the Kruskal-Wallis test and if significant, then multiple Wilcoxon-rank test. We also examined shape sexual dimorphism using Principal Component Analysis, Thin-Plate-Spline, and a MANOVA test, as described previously.
Multivariate regressions of shape variables on log-transformed centroid sizes were used to analyse allometry.

Size variation considering the four subspecies altogether
Forelimb proximal phalanx. Kruskal-Wallis tests show significant size differences in almost all of the variables: subspecies (p<0.002), habitat (p<0.05), mobility (p<0.05), and most of the interactions between these variables (Table 6). Mobility+Habitat interaction did not yield significant differences. When subspecies category is included in the interaction for both mobility and habitat, there are always significant differences (p<0.05). However, it must be noted that, for all subspecies expect R.t. caribou, only one type of habitat (tundra) and mobility (planarly) are represented in our sample, a composition that affects the statistical Pairwise comparisons did not show size differences when they were displayed by subspecies +mobility+habitat. However, although it is not statistically significant, larger sizes correspond to male specimens which performs altitudinal migrations (R. t. caribou, mountain setting) ( Table 7) as it can be observed in boxplots (Fig 3).

Size variation within Rangifer tarandus caribou
R. t. caribou is the only subspecies in our sample that include specimens of varying habitats and mobility. Kruskal-Wallis yielded no statistical differences (Table 9) on size for Rangifer tarandus caribou subspecies. Therefore, no pairwise comparisons were carried out. A very small size variation is perceptible in the boxplots (see S2 File) when both habitat and mobility are considered, even if specimens that carried out altitudinal migration have slightly bigger phalanges. These small differences could be explained, at least in part, by differences in sex ratios between subpopulations (see Table 2), although this would be difficult to explore without a larger sample size.

Shape variation results considering the four subspecies altogether
Forelimb proximal phalanx. The first 2 PCs explain 46.42% of the morphological variance ( Fig 5). Shape variation along PC1 accounted for 28.65% of the variance while PC2, 17.77%. Although there are overlaps according to both Habitat and Mobility type (Sedentary, habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal) using Wilcoxon rank sum test after "Benjamini-Hochberg" correction. A significant contribution was considered for P value <0.05 (in bold). https://doi.org/10.1371/journal.pone.0285487.t007 Table 8. Pairwise comparisons of log-transformed centroid-sizes. Pairwise comparisons of log-transformed centroid sized for hindlimb phalanx between the different groups according to subspecies Towards the positive values of PC1, proximal morphology is more massive towards cranial side and mediolaterally narrower (Fig 6). Axial ridge resulted to be straighter and abaxial articular surface broader. However, towards negative values, proximal joint appeared to be more prominent mediolaterally than craniocaudally. Middle groove was deeper towards the negative values.  Thus, the abaxial view is more massive towards the positive values of PC1 (mainly individuals from tundra habitat and planarly migration) and narrower towards the negative values of the same axis (mostly Rtcaribou_Mountain_Altitudinal specimens). In the caudal view, the morphology tends to be more inwardly curved towards positive values of PC1 and straighter towards negative ones. In the same way, in the axial view, a greater protuberance is observed in the lower part of the articular surfaces towards the positive values of the PC1 axis and a flatter configuration towards the negative ones.

Subspecies
Regarding PC2, the main differences are observed in the proximal view of the articular surface. In the axial view, a narrower configuration is observed towards the positive values of this

PLOS ONE
axis. The theoretical shape at the PC2 minimum showed a proximal joint that is larger craniocaudally and narrower on an axial-abaxial axe and the opposite it is shown towards its maximum shape and occurs more often on mountain specimens.
CVA results for classification (Table 12) obtained 82.35% overall accuracy according to habitat and a lower percentage in accuracy when it was performed according to mobility pattern (73.52%). These percentages are a bit higher when calculated according to 2 mobility groups (migratory/sedentary), 82.3%. Thus, classification accuracy in this case improves when mobility patterns are grouped into two categories rather than multiple categories. In all cases, Kappa statistic (see S1 File) was higher than 0.6, which indicates a substantial agreement beyond chance.

Hindlimb phalanges
First two PCs accounts for 48.61% of the morphological variance (Fig 7). Rtcaribou_BFor-  The PC1 minimum configuration (Fig 8) seemed to be more stretched craniodorsally and narrower mediolaterally. Contrary, towards the positive values of PC1 or maximum configuration expressed a theoretical morphology that was broader mediolaterally and more elongated dorsally in the lateral articular surface. Towards positive values of PC1 (Fig 8), the morphology tends to be more curved towards the abaxial side. On this same side, the configuration is more massive towards positive than negative values, however in the caudal and axial views the morphology is wider and more pronounced towards negative than positive values of the PC1 axis.
In PC2 (Fig 8), the morphological configuration in this axis stands out for being more pronounced towards its positive values. The abaxial view shows a greater protuberance in the lower part of the articular surfaces towards the positive values of the PC2 axis and a flatter configuration towards the negative ones. However, in axial view, it tends to bulge more towards negative values, while being flatter towards positive values of the axis.

PLOS ONE
Reconstructing past migratory behaviour of reindeer (Rangifer tarandus) Regarding PC2 (Fig 8), we find the opposite as in PC1. The theoretical morphology at the minimum appeared to be more massive mediolaterally in the proximal articular surface and with a deeper groove. However, it was narrower abaxially at the PC2 maximum.

PLOS ONE
Reconstructing past migratory behaviour of reindeer (Rangifer tarandus)

Shape variation amongst Rangifer tarandus caribou
In forelimb phalanges, the first two PCs account for 49.76% of the total variance ( Fig 9A). Individuals with positive PC2 values (sedentary specimens) display a joint morphology (Fig 10) that is more stretched abaxially and elongated in the craniocaudal axis; however, individuals with negative PC2 values have a joint that is wider in the abaxial axis (Rtcaribou_Tun-dra_Planarly and Rtcaribou_Mountain_Altitudinal).
In the case of hindlimb phalanges, the analysis of the PCA reveals a complex configuration. PC1 and PC2 account for 50% of the variance (Fig 9B). Rtcaribou_tundra and Mountain are mainly distributed along PC1 with positive and negative values. In contrast, specimens of Rtcaribou_BForest_Sedentary are mainly distributed along positive values of PC1 and PC2, overlapping with some Rtcaribou_Mountain_Altitudinal individuals.
Thus, specimens toward positive values of PC1 are shown as having an axial articular surface that is more elongated than those towards negative values (Fig 11). Those specimens towards PC2 positive values show a stretched morphology on the abaxial axis and a more elongated axial surface towards the caudal side, whereas those that have negative values have the opposite configuration (wider abaxially and narrower craniocaudally).
In both fore-and hindlimb phalanges (Table 14), there are significant differences (p<0.05) according to habitat, mobility, and the interactions between the two. Nonetheless, pairwise comparisons revealed differences between all habitat categories and between altitudinal and planarly migration types in forelimb phalanx (Tables 15 and 16).

Sexual dimorphism on size (R.t. caribou)
Despite the fact that sex is frequently considered an important variable in reindeer morphological studies [75,77,83,88] size and shape sexual dimorphism can hardly be addressed in our

PLOS ONE
Reconstructing past migratory behaviour of reindeer (Rangifer tarandus) sample, due to the low number of specimens of known sex ( Table 2). Only the R. t. caribou sample have specimens of both sexes. The lack of enough specimens prevents further analyses due to low statistical power to detect meaningful differences between groups. In contrast, boxplots were used to visualize log-centroid size distribution. Thus, when the variation in logtransformed centroid size is displayed in boxplots (see S2 File) according to male/female it can be observed that males tend to be larger in size than females, although there is an important overlap. Considering these results and our limited sample in terms of specimens of known sex, we discarded the sex variable for further analyses.

Allometry
For forelimb phalanx, the percentage of variance related to size was low, only a 5.90% and the allometry was significant (p = 0.04). However, for hindlimb phalanx, the percentage of variance related to shape was even lower, at 2.34%, and allometry was not significant (p = 0.38).
Multivariate regression (Fig 12) of shape scores against log-transformed centroid size showed a clear separation of R. t. caribou with the rest of the subspecies (except R. t. granti) in forelimb phalanx. For hindlimb phalanx (Fig 12), R. t. peary is completely separated from the rest of the subspecies, which overlap among them.

Discussion
In the present study, we have examined the morphology of caribou proximal phalanx and its relationship with mobility and habitat type in order to potentially apply into archaeological

PLOS ONE
record for palaeoecological reconstructions. Despite the fact that a wider and more complete sample for further analyses is needed, this work represents a first step in this regard.
As we have already mentioned along our study, the main limitations of our sample are the size and the lack of information on the sex of the specimens. Due to the fact that sedentary populations are classified as "Endangered", the acquisition of samples from these herds was more difficult (since they are barely extant in museums, it was only possible to acquire them through samples sent by biologists acquired when a caribou died in a traffic accident or was

PLOS ONE
Reconstructing past migratory behaviour of reindeer (Rangifer tarandus) found dead in the forest), resulting in uneven sample sizes. Regarding sex, most of the specimens from museum collections were collected at the end of the 19th or beginning of the 20th century and they did not specify the sex, so we can only know it from those that were collected in situ by biologists nowadays. Because only few individuals within the sample population are of known sex, we cannot directly evaluate the relationship between tested values and sex regarding to habitat or mobility. Specimens with bone pathologies [71,72] were excluded from the study to avoid this bias in our sample. Despite the limitations imposed by our sample size, there is a consonance between the statistical analyses and the results in the configuration of the morphospace. Due to this reason, we considered including the subspecies category as a preliminary step while waiting for a larger sample size.
Studying joint surfaces allows for the use of complete or fragmented bones, increasing the sample size that can be used in archaeological records [67]. Furthermore, when compared to previous studies with linear measurements on reindeer phalanges, a 3D GMM study has the advantage of allowing a deeper study in the morphology of these phalanges, and is able to accurately capture shape variation, even with a small sample size (e.g. altitudinal migration was difficult to explore using linear measurements, cf. [41]). For validation of the methods presented here, we recommend similar studies on caribou populations in other regions so that we can assess the preliminary variations we observed.

Interpreting size variation: A complex interplay of sexual dimorphism and migration costs?
Size variations are mainly imputable to subspecies differences. R. t. caribou are statistically significantly larger than R. t. pearyi, which was expected, since the latter is amongst the smallest R. t. subspecies [119]. Peary Caribou, inhabit the Canadian Arctic Archipelago (island syndrome) and northern of the Boothia Peninsula [94,119]. They are adapted to the Arctic environment having a compact body size for conserving heat (Allen´s rule) [120] and to the limited plant growth with a highly compressed growing season and long periods of snow-covered frozen standing vegetation [119]. All these facts explain their lower centroid sizes.
At an intra-subspecies scale, size differences are not significant between populations of R.t. caribou. Forelimb phalanges of specimens that perform altitudinal migration tend to be slightly larger than those belonging to planarly or sedentary migratory behaviour, but it is yet unknown if this is linked to real population differences or unbalanced sex ratios in our sample.
If confirmed, these size differences might be related to the migration costs. Migration of long distances might increase access to high quality food or reduce the risk of predation but it also involves an energetic cost that may be reflected in body conditions of animals [121,122]. The proportional cost of horizontal or planarly travel that tundra caribou performs increases with decreasing body size and long movements, particularly in winter, have a negative impact on caribou body size [123,124]. However, at present, there is no evidence of a statistically significant influence of habitat or mobility on phalanx size.
Size differences tend to be more pronounced on forelimb than hindlimb members. This could be explained by the fact that forelimb supports more stress as they carry greater share of the body weight [75][76][77][125][126][127][128]. The use of forelimb more actively is also linked to the search of food, breaking through ice layers in many cases [69,70,75]. It is, however, important to keep in mind that these results do not satisfy the thresholds of statistical significance, and therefore should be interpreted carefully.

Shape variation: The impact of subspecies, substrate type and snow cover?
Despite the fact that significant shape differences exist between populations, it is difficult to identify a specific morphology that corresponds to a given habitat or mobility type.
If we consider CVA, in terms of classification instead of associated morphological description, CVA is specifically designed to identify linear combinations of variables that best discriminate between groups. It maximizes the separation between groups while minimizing the variation within groups, and thus can be more effective in identifying the underlying differences between groups. For forelimb phalanges, the results suggest that habitat may be a stronger predictor of classification accuracy than mobility pattern. However, on hindlimb phalanges, these ones indicate that both habitat type and mobility pattern are important predictors of the observed patterns. Thus, these outcomes suggest that the morphology of caribou forelimb phalanges is more heavily influenced by the type of habitat in which they inhabit than their mobility pattern. This could mean that the forelimbs are adapted for specific types of terrain (e.g. rocky and mountainous terrain). In contrast, the results for the hindlimb phalanges indicate that both habitat type and mobility pattern are significant predictors of the observed patterns, suggesting that they are better suited for more dynamic activities (e.g. running). Thus, the different predictors of classification accuracy between the two limb sets likely reflect the different functional roles of each set of phalanges. The results are consistent with the MANOVA results, which also indicate that habitat variable plays an important role in the configuration of the shape of the phalanx morphology.
Results on forelimb phalanx revealed, along PC1, on one hand, a joint with a more massive morphology towards cranial side and medio-laterally narrower (mainly present in Rtcari-bou_BForest_Sed and partially on Rtcaribou_Tundra_Plan), and on the other hand, proximal joints that appeared to be more pronounced mediolaterally rather than craniocaudally (mostly Table 19. MANOVA pairwise comparisons results (considering all the categories). There is however no clear distinction according to habitat or mobility. In spite of the difficulty of establishing the functional implications of these morphological differences, it could suggest adaptations to different types of movement or stresses placed on the forelimb phalanges in different environments. For example, the more massive joint towards the cranial side and medio-laterally narrower morphology observed in individuals inhabiting the boreal forest may provide greater stability when moving through snow, whereas the more pronounced mediolateral morphology observed in mountain individuals may provide greater flexibility when traversing uneven terrain. Towards the maximum configuration of PC2 we observe a morphology that is wider in the axial-abaxial axe and narrower craniocaudally (especially for sedentary R.t. caribou from the mountains, and also for planarly migrating R.t. granti from the tundra) and the opposite configuration along negative values of PC2, with specimens that are both sedentary or migratory, inhabiting mountain, boreal forest or tundra. At this point, PC2 is playing a smaller role in determining the shape of the joint., In other words, there is no clear distinction between the morphology of populations based on their migration behaviour or habitat.

Rtcaribou_BorealForest_Sedentary-Rtcaribou_Mountain_Altitudinal
Forelimb phalanx from sedentary R.t. caribou inhabiting the mountains overlaps to a lesser extent in the morphospace with the other populations (especially when the R.t. caribou subspecies is analyzed separately) and reflects a more characteristic morphology, contrasting with the boreal forest population, which is also sedentary. The mountain caribou's forelimb phalanx is likely to be more robust and adapted to their mountainous environment, whereas the boreal forest population's forelimb phalanx is likely to be more adapted to the conditions found in the boreal forest (e.g. a more varied terrain, with greater amounts of wetland; a more diverse array of vegetation, including shrubs, grasses, and forbs; and a greater variety of food source). In both cases, it is worth clarifying that when we refer to sedentary populations we do so to talk about those that travel distances less than 200 km annually [see 41], so in this case substrate and snow coverage potentially have an important impact on the morphological configuration of the joint of the proximal forelimb phalanx. Our sample of boreal forest caribou comes mainly from Québec (Canada), where the forests are covered in snow more than half of the year [88]. These types of boreal forest populations also seek places with snow that can be used as refuges from predators. Sedentary montane populations, however, remain in alpine areas to feed on lichens in windswept areas where snow is absent or shallow [48]. These differences in substrate type and snow cover may also explain the slight overlap in morphological configuration between sedentary and altitudinally migrating mountain individuals: even though both inhabit mountainous ecosystems, they move in different ways (altitudinal migration vs. sedentary). While for the sedentary mountain caribou we have already discussed the characteristics of their habitat, for those that perform altitudinal migration, at the end of autumn, when the snow begins, they migrate from upper to lower elevations where the snow is shallow enough, avoiding open habitats where the snow is harder, again beginning its travel to higher elevations around February towards subalpine forests [48,129]. This suggests that the functional morphology of this particular subspecies (Rtcaribou_Mountain_Sedentary) is well adapted to its environment, allowing to thrive in its mountainous habitat (characterized by a diverse range of elevations, and the presence of high-altitude alpine tundra, ranging from rolling hills and valleys, to steep rocky slopes, and high peaks).
Morphological differences are even harder to classify on hindlimb phalanges. All populations are distributed along PC1, and PC2 poorly distinguishes populations, indicating that the morphology of the hindlimb phalanx of each subspecies is largely dependent on the habitat and mobility pattern of the subspecies. As a result, the morphological variation of hindlimb phalanx is not so evident to interpret.
MANOVA tests showed statistically significant differences in relation to Habitat for both fore-and hindlimb phalanges, and in relation to Subspecies and Mobility for hindlimb phalanx. Shape differences were more pronounced on hindlimb than forelimb phalanges, and more marked when pooled by Habitat in both cases. This may suggest that the habitat in which Rangifer tarandus lives, has a significant impact on the joint shape of forelimb first phalanges, as has been proven in different Bovidae species [60, 67,130]. To understand habitat´s influence, it is important to note the variety of environments and substrates that caribou inhabit, ranging from closed boreal forest to treeless arctic tundra passing through the mountain with a mosaic of landscapes ranging from open alpine areas to subalpine areas and lowelevation forests [69]. There is no doubt that the snow adaptation impacts the shape of the proximal joint, not only during migration, but also when the front hooves are used to dig through the snow to find lichen during snowy seasons [68].

Shape variability at the intra-subspecies scale: The key influence of habitat?
Based on pairwise comparisons, in hindlimb phalanx, main differences on shape are found between the different subspecies, a fact that motivated the realization of the same analyses within the Rangifer tarandus caribou subspecies. Such intra-subspecies analysis confirmed that phalanx shape is statistically different according to Habitat and Mobility, even if only R.t. caribou specimens are considered.
When R. t. caribou was examined separately, MANOVA showed differences in habitat and mobility, as well as their interaction, in fore-and hindlimb phalanges. In pairwise comparisons, all habitat categories showed significant differences, regardless of whether they were fore-or hindleg and demonstrating (as determined by the p values) how crucial the impact of habitat on shape is. Regarding mobility, there were differences for forelimb phalanges between planarly on one hand, altitudinal and sedentary on the other, and only between altitudinal and planarly in the hindlimb ones. Considering the interaction, pairwise comparisons did not yield significantly different results in foreleg phalanx, but in hindleg, between Rtcaribou_BFor-est_Sed and Mountain_Altitudinal and Tundra_Planarly (respectively). Consequently, habitat has again the greatest impact on the shape, although to a lesser extent the type of mobility also has an impact. The altitudinal migration in the mountains occurs vertically, through rugged and mountainous substrates in response to snowfall conditions, forage availability, and to avoid predators. Migration in tundra occurs horizontally, with a characteristic vegetation of dwarf and upright shrubs [94,129,131]. These differences are not only in the substrate and the type of movement, but also in the distance traveled by the herds. Rtcaribou_Tundra_Planarly travels distances between 1,220 and 1,770 km annually, while Rtcaribou_Mountain_Altitudinal can travel between 140-240 km only in the spring migration [81,90,92,131]. Furthermore, the marked difference in the hindlimb phalanx may also be explained by the implications of propulsion (along with what was previously discussed), resulting in a higher level of stress among individuals who move vertically instead of horizontally due to the substrate conditions, in addition to the fact that in quadrupeds, as it has been seen in other species, e.g. rhino, horse [126,127,132] the hind limbs bear relatively less weight than the fore limbs and play an additional propulsive role during locomotion.
The results show that subspecies structure the data significantly, especially for size, but habitat and mobility also have significant effects on phalanx morphology. Shape analyses on Rangifer tarandus caribou subspecies yielded shape differences that are especially significant in both fore and hindlimb phalanges according to habitat or mobility categories. Consequently, if the subspecies variable is excluded, we can examine how habitat has a greater impact on the morphological configuration. This can be explained by the wide variety and different habitats used by Rangifer tarandus caribou: those that inhabit mountain environments move among windswept alpine ridges, subalpine forests, and low-elevation forests [69]. In the Boreal Forest, in which the snow is a key feature, they seek habitats with thin and soft snow cover with a more sedentary lifestyle [88]. Those inhabiting the tundra will find many types of flowering plants, mosses, and lichens, as well as shrubs, grasses, and sedges which form swards and tussocks scattered across the tundra [94]. Thus, e.g. Rtcaribou_BForest_Sed tends to have a forelimb phalanx joint morphology with a more cranially prominent axial articular surface, compared with sedentary mountain individuals. This morphology could be linked to the adaptations that enable it to better handle the specific demands of the boreal forest environment (e.g. deep snow and the ability to dig through the snow to reach vegetation; rugged terrain and extreme cold). For example, the cranially prominent axial articular surface may allow for a greater range of motion in the joint, or it may provide additional surface area for muscle attachment, allowing for more precise and powerful movements. However, further research would be necessary to fully understand the functional implications of this morphological difference.

Preliminary distribution on sexual dimorphism
Even though the sample size by sex is limited for both forelimb and hindlimb phalanges, there are no differences between male and female specimens from the same subspecies (R. t. caribou) when log-centroid sizes were displayed on boxplots to explore their distribution. In this regard, Hull [71] noted in her study to distinguish forelimb and hindlimb phalanges, that, even though differences between sexes were appreciated, they also included a large overlap, so she concluded they cannot be included in the study without additional context. The boxplots show that male phalanges are slightly larger than female ones, but they overlap. These non-significant results do not imply that these differences between the two sexes do not exist in our sample, as has been demonstrated for other anatomical parts [82], it can simply reflect that phalanges are not an indicator of this characteristic and/or a larger sample and further studies are needed.

Allometry
Allometry is also more significant on forelimb than on hindlimb (absent in our sample) phalanges, similar to what has already been observed in other works. Previous studies on reindeer [75,77] show that allometry was generally stronger on the forelimbs than the hindlimbs as well. It has been observed in other species as well, for example the rhinoceros [127,128] and some bovids [130]. The fore-phalanx allometry is not very significant, but it may be compensated by the rest of the forelimbs. Thus, a predictable morphological and biomechanical response of the bone is necessary to resist the stresses caused by body mass, which is heavier (as noted above) on forelimbs due to antlers (which are heavier in males), resulting in additional forces and constraints on the foreleg bones.

Conclusion
Reconstructing the interactions between Palaeolithic human groups and reindeer is often complicated by a general lack of data on the ethology and ecology of these past animal communities. At the end of the Late Glacial, reindeer gradually disappeared from southwestern Europe with the climate warming. As they became more fragmented geographically, different responses should be expected, such as variation in their migratory patterns. As the migration behaviours of reindeer herds during this period are largely unknown, it is difficult to interpret Late Glacial socio-economic organization and mobility strategies.
Due to the fact that locomotor morphological adaptations are closely related to habitat preferences and mobility [67,78,133,134], it is possible to infer behavioural responses using an animal's cranial and post-cranial morphology, in the case of our study, using first phalanx. Besides, phalanges are one of the most useful bones for palaeocological purposes because these elements are commonly preserved in fossil assemblages. Together with 3D geometric morphometrics, we were able to examine their morphology from an anatomical region of interest and visualize morphological variation. The present study has contributed to our understanding of the caribou proximal phalanx's functional morphology and its relationship to mobility and habitat type (adapted to provide support on uneven terrain and harsh habitats), but it still cannot be easily applied to archaeological bone remains. A complex interplay of factors has been shown, with the impact on phalanx morphology of subspecies status, sexual dimorphism, substrate type, snow cover, migration costs, habitat, etc. This diversity of interrelated factors makes interpretation difficult with the limited sample of modern reindeer phalanges currently available. Continued efforts in the 3D scanning of phalanges should lead to further progress.
Although future studies with larger sample size are needed, this contribution represents a necessary first step to explore reindeer phalanx morphology and size differences for potential archaeological applications.
Supporting information S1